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Abstract 



ON 



Multistep denaturation in a short circular DNA molecule is analyzed by a mesoscopic Hamilto 
nian model which accounts for the helicoidal geometry. Computation of melting profiles by the path 
integral method suggests that stacking anharmonicity stabilizes the double helix against thermal 
disruption of the hydrogen bonds. Twisting is essential in the model to capture the importance 

of nonlinear effects on the thermodynamical properties. In a ladder model with zero twist, anhar- 

g 

monic stacking scarcely affects the thermodynamics. Moderately untwisted helices, with respect 

to the equilibrium conformation, show an energetic advantage against the overtwisted ones. Ac- 
cordingly moderately untwisted helices better sustain local fluctuational openings and make more 

i-H unlikely the thermally driven complete strand separation. 
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1. Introduction 

In the transcription process, the RNA polymerase attaches to the DNA near the start of 
the gene and opens a small segment of the double helix thus reading and coping the coded 
instructions into a template, the mRNA, for the synthesis of proteins. Both the molecular 
machine and the transcription bubble move along the gene by activation processes associ- 
ated to proteins which bind to DNA portions, enhancers p] , located even several thousands 
bases away along the chain with respect to the promoter starting region [2]. In some cases, 
long range enhancer-promoter interactions appear to be mediated by DNA structural defor- 
mations such as looping, bending and twisting whose molecular origin is under investigation. 
Certainly the flexibility of the double helix favors long range effects in biological functioning 
[3H5] whereas twist deformations modify the three dimensional structure, producing com- 
pact supercoiled shapes which are crucial for DNA packaging within the nucleus [6H9]. It is 
known that the opening of transcription bubbles is possible only if accompanied by a local 
decrease in the twist angle with respect to the equilibrium value. Unwinding is measured 
by a negative supercoiling and the undertwisting of DNA can yield concentration of energy 
such to break hydrogen bonds, eventually leading to denaturation [TQl UT] . Twisting is de- 
scribed by the angle that the base pairs rotate around the molecule axis. B-DNA at room 
temperature has a helix repeat of ~ 35 A hosting ~ 10.4 base pairs, hence the equilibrium 
twist angle is 9 eq ~ 0.6 rad. As proposed by Prohofsky [12], such helical structure may be the 
optimal configuration to capture the chemical energy flux which nourishes RNA polymerase 
in DNA transcription. Similar arguments have been recently raised by a numerical study of 
the nonlinear dynamics in helicoidal DNA showing that the energy localization required for 
bubble formation is driven by anharmonic stiffness [13]. In this regard, a microscopical un- 
derstanding of 9 eq in unperturbed DNA seems therefore preliminary to the comprehension of 
those local fluctuational openings occurring in replication and transcription events [T4"l [To] . 
As the latter involve disruption of the base pairs hydrogen bonds in a way similar to the 
denaturation process, whether mechanically or thermally driven [16HT9], the properties of 
the melting transition have drawn such a great interest over the last decades [20Tl3"l~] . 

To elucidate the DNA biological functions in terms of the fundamental interactions, path 
integral techniques [32l [33] are applied in this paper to the thermodynamics of a circular 
molecule described by a mesoscopic nonlinear Hamiltonian, the Dauxois-Peyrard-Bishop 



(DPB) model [3U [35] . Hamiltonian approaches have the advantage of treating the system 
at the level of the base pair [SB] thus permitting to study the thermodynamic behavior by 
tuning the heterogeneity degree in the sequence. 

In a previous investigation [37], the standard DPB model has been enriched by intro- 
ducing a rotational term in the stacking potential thus accounting for the twist between 
adjacent base pairs along the molecule axis. It has been found that the equilibrium twist 6 eq 
provides in fact the stablest conformation against thermal disruption of the base pair bonds 
in the denaturation regime. As the twist on its own gives stability to the system, it has con- 
sistently been shown that the melting profiles of the twisted conformation are less affected 
by changes in the harmonic backbone coupling than the untwisted DPB model. However 
no analysis of the anharmonic stacking effects has been carried out in Ref.[37j. As the 
DPB model attributes a special role to the anharmonic backbone interactions in sharpening 
the denaturation transition [38], here I study specifically the relation between anharmonic 
stacking and twist angle in a modified DPB model which also includes a stabilizing solvent 
interaction term. It is found that the twisting of the double helix is essential to switch on 
the anharmonic stacking effects on the melting profiles whereas stacking anharmonicity is 
negligible in the standard DPB model. 

The system Hamiltonian is presented in Section 2 together with a discussion of the stack- 
ing potential energetics as a function of the twist. Section 3 outlines the path integral 
computational method. Section 4 contains the results: the melting profiles are displayed 
both for underwound and overwound short double helices while the effects of stacking an- 
harmonicity are reported for the equilibrium twist case and for the untwisted DPB model. 
Also the role of the solvent potential is discussed. Some final remarks are made in Section 
5. 

2. Model 

To model a double helix of iV base pairs with reduced mass //, I start from a generalized 
DPB Hamiltonian in which the pair mates separation, with respect to the ground state 
position, is given by the transverse stretching y n (for the n-th base pair). These are generally 
much larger than the longitudinal base pairs displacements along the molecule backbone 
which are accordingly dropped. Thus, the one-dimensional Hamiltonian reads: 
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V M {Vn) = D n (exp(-a n y n ) - l) 
K 

V S (y n , Vn-l) = —g n ,n-l(yl ~ ^VnVn-1 COS 6» + yl_ x ) 

g n , n -i = 1 + pexp[-a(y n + y n -i)] 

V S oi(y n ) = -D n f s (t&nh(y n /l s ) - l) . (1) 

Like the DPB Hamiltonian, Eq. (IIJ incorporates nonlinearities both in the hydrogen bond 
base pair interactions shaped by the Morse potential Vjy(y n ) and in the stacking coupling 
Vs(y n ,y n -i) between neighboring bases along the two strands. D n and a n are the pair 
dissociation energy and the inverse length setting the hydrogen bond potential range: they 
account for heterogeneity in the sequence. 

Considering that: i) the AT- and GC-base pairs are linked by two and three hydrogen 
bonds respectively, ii) there is electrostatic repulsion between the highly charged phosphate 
groups on complementary strands, the pair dissociation energies can be reasonably taken 
as Dat = SOmeV and Dec = 4:5meV. These values are above fc^T at room temperature. 
ks is the Boltzmann constant and T is the temperature. Larger energy values [39l HO] 
would shift the denaturation steps at higher T. The inverse lengths are a^r = 4.2A and 
a cc — 5A in view of the fact that AT-base pairs have larger displacements than GC-base 
pairs. 

K is the harmonic stacking related to \x by K = pu 2 , v being the frequency of the 
phonon mode, p and a are the anharmonic stacking parameters which are also assumed 
independent of the type of base at the n and n — l sites. The homogeneity assumption [H] 
for the stacking potential relies on the observation that both types of base pairs contain 
a purine plus a pyrimidine, the former being larger and heavier. Thus the AT- and GC- 
base pairs are comparable in size and weight with their effective masses usually taken as 
\x ~ 300a. m. u.. While the inclusion of stacking heterogeneity is motivated in a growing 
number of studies [121 - H5] and may be straightforwardly incorporated in the path integral 
approach, there is no need to further enlarge the set of input parameters to the purposes of 
the present work. 
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I take values in the commonly used range, K = 60meVA , p = 2, a = 0.5A whereas 



the anharmonicity effects due to larger p are discussed below. Somewhat smaller (or larger) 
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harmonic couplings also found in the literature, i.e.: K ~ 25meVA [39J, fT ~ lOOmeVA 



would have the main effect in our method to shift downwards (or upwards) the denat- 
uration temperature while leaving unaltered character of the transition and shape of the 
melting profiles [47] . 

Beyond the DPB Hamiltonian, Eq. (fl| introduces: 1) Torsional effects [37J EH] through 
the angle 9 between two adjacent base pairs at the n and n — 1 sites along the molecule 
backbone. While the stacking is even function of 6, only positive twists are allowed to 
describe right-handed helicoidal structures and to avoid zig-zag configurations which may 
also minimize the energy. 2) The solvent potential V so i(y n ), as proposed in Ref. [32J 150] . 
which adds to the Morse term thus enhancing by f s D n the height of the energy barrier 
above which the base pair dissociates. 

The solvent factor f s mimics the effect of the salt concentration [iVa + ] on the Morse 
potential depth D n . The latter is known to scale essentially linearly with the melting tem- 
perature T m of homogeneous GC or AT sequences (whereas, in the harmonic strong coupling 
limit, T m oc \/D n ). As T m varies logarithmically with [A^a + ], also the dissociation energy is 
assumed to increase as D„ oc ln[Na + ] [51]. Accordingly, I take f s = 0.3 simulating a high 
counterion concentrations in the solvent, [A^a + ] ~ 1M , which screens the negatively charged 
phosphate groups j52j|53]. The length l s defines the range beyond which the Morse potential 
plateau is recovered and D n returns to be the fundamental energy scale. For y n > l s , the two 
strands are apart from each other and the hydrogen bond with the solvent is established. 
Hereafter I take l s = 3A whereas the possible effects due to the extension of the solvent 
potential range will be discussed in Section 4.3. 

While Eq. (fTl) refers to a system with N + 1 base pairs the presence of an extra base pair 
2/o is usually handled by taking periodic boundary conditions, y = yjy, which close the finite 
chain into a loop. This makes the model suitable for application to circular DNA although 
bending effects [HU [55] associated with tilting between base pair planes are not considered 
in this paper. 



2.1 Stacking versus Twist 

I propose an intuitive argument, based on the principle of the energetic convenience, to 
support the idea that a relaxed DNA molecule modeled by Eq. ([l]) should plausibly prefer 
twist configurations with 9 ~ 9 eq . To begin with, let's assume that: i) N = 100 and ii) 
in the equilibrium case, the number of base pairs per helix turn is h = 10. As N/h is an 
integer and the helix axis lies in a plane, the strand ends line up precisely so that there is 
no need to introduce a curvature which would generate some writhe [56H58] . In this case, 
the linking number coincides with the twist, Lk = Tw, the equilibrium twist of the relaxed 
configuration is (Tw) eq = N/h = 10 and 9 eq = 2n/h. 

In replication and transcription, DNA unwinds thus partly dissipating the torsional stress 
associated with the supercoiled geometry Double helix unwinding amounts to a twist reduc- 
tion with respect to (Tw) eq and generates negative supercoiling which is likely an efficient 
strategy to fluctuate in solution. While supercoiling generally refers to the molecule axis 
coiling upon itself, in our simple model the term indicates a variation in the helical repeat. 
Thus the model in Eq. (IT]) should also account for an energetic advantage of underwound 
configurations, Tw < 10, compared to the overwound ones, Tw > 10. To tune the twist, 
imagine to cut the molecule, stretch the double helix and change the helical repeat so that 
Tw is an integer before re-closing the ring. The DPB ladder model corresponds to the limit 
case, h — > oo [55] . 

The intra-strand and inter-strand interactions, described by Vs(y n , y n -i) and VM(y n ) 
respectively, compete on the energy scale. If, for some displacements y n , either one or the 
other (or both) potential term is very large the corresponding contribution to the partition 
function becomes vanishingly small. Vm^m) is in principle unbound on the y n < semiaxis 
where negative displacements depict the two complementary strands coming closer to each 
other than in the equilibrium state, y n ~ 0. The negatively charged phosphate groups 
prevent however the two strands from coming too close and the hard core of Vm (y n ) accounts 
for such repulsive effect. Thus the base pair transverse stretchings encounter a physical cutoff 
y m in which is set in our computational method [33] by the requirement VM(y m in) ~ D n . For 
AT-base pairs, it is found that y m in ~ —0.2A while y n <C y m in would lead to inconveniently 
large electrostatic repulsions. Likewise, one selects the appropriate range for GC-base pairs. 
Instead, on the y n > semiaxis Vuiyn) is bound hence, arbitrarily large displacements may 
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FIG. 1: (Color online) Stacking potential as a function of the difference between the displacements 
of two adjacent base pairs along the molecule backbone. Four twist conformations are assumed: 
(a) Dauxois-Peyrard-Bishop ladder model; (b) Strongly underwound; (c) Equilibrium twist; (d) 
Strongly overwound, p is dimensionless; a is in units A ; K is in meVK~ l . 



be included in the computation. Nonetheless, once the base pair dissociates, the pair mates 
tie to the solvent thus there is no reason to further extend the range of the stretchings. 
Moreover, very large y n would yield constant energy contributions to the partition function 
with no effect on the free energy derivatives. A cutoff y ma x ~ h is then generally appropriate 
while, to our present purposes, even y max ~ 0.8A suffices. 

With these premises, I discuss in Fig. [l] the stacking potential energetics as a function 
of the twist conformation tuned by the helical repeat. Two adjacent base pairs along the 
molecule backbone, y\ and 7/2, are assumed to vary in the range [—0.2, 0.8]A and the relative 
stacking Vs{y2,yi) is plotted versus the difference y± — y 2 pinning y x to six values in that 
range. 

Fig. Hi a) refers to the DPB ladder model, 6 = in Eq. pj). Whatever value is assigned 



to 2/1, the second base pair can adjust its position so as to keep Vs low in energy: i.e., setting 
Hi = 0.8 A and varying y 2 , one finds a range y\ — y 2 G [0, 0.5] A in which Vs does not exceed 
~ 30meV, the Morse plateau for AT-pairs. Assuming y\ = OAA, Vs < 30meV in the broad 
range y\ — y 2 G [— 0.4, 0.4]A. Likewise for any other y\. Obviously the potential minimum 
always occurs for y 2 = y\ but the stacking energies remain generally of order of Dat even 
in the case of large differences y± — y 2 . Because of a fluctuation, a base pair may get a 
significant displacement (~ 0.8A) out of the equilibrium nevertheless the adjacent base pair 
can move quite freely keeping the stacking potential at values comparable with Vm- Thus the 
stacking of the DPB model gives little stability to the double strand molecule as it does not 
sufficiently discourage large relative fluctuations which therefore contribute to the partition 
function. As a general requirement the model should be both flexible enough to absorb 
fluctuations and stacked enough to guarantee some stability against thermal disruption. 

This goal is achieved by introducing a twist through a finite 9 as shown in Fig. [lib) 
and, thoroughly, in Fig. [11(c) which displays the equilibrium conformation. Mostly in the 
latter, according to the value assigned to y%, some y 2 displacements are discouraged as they 
appear too costly on the energy scale. If y\ ~ 0.8A, y 2 can only assume large enough values 
(such that y\ — y 2 G [0,0. 3]A) to minimize the stacking energy but, even in this restricted 
range, Vs is much larger than in Fig. pTa) signalling that the torsion is indeed effective in 
providing some restoring force to the molecule. The energy spacings among the minima in 
Fig. [Tj(c) indicate that some couples of base pairs (j/i, y 2 ) are less likely to contribute to the 
partition function and, at the same time, they also show the system capability in accepting 
moderately large fluctuations. Taking y\ = OAA, the possible displacements for y 2 which 
cost less than 30meV are such that y\ — y 2 G [~ — 0.35,0. 4] A hence, y 2 has a broad range 
to explore and can become even larger than y\. Unlike Fig. HTa), the potential minimum 
does not correspond to the translation y 2 = yi but rather to a slightly positive y\ — y 2 . y 2 
follows y\ to pull it back. This is the restoring effect yielding flexibility to the system. 

One may wonder whether larger twists, overwound conformations, might become even 
more energetically convenient than Tw = 10. The answer is however negative as shown in 
Fig. [lid) where Vs is plotted for the case Tw = 15. Now the energy scale changes drastically 
and the energy cost to follow a sizeable fluctuation in y\ becomes too high: if y\ > OAA, 
there is no y 2 value to keep Vs comparable to the Morse dissociation energy. Only taking 
yi in the range [—0.2, 0.2]A there exists a likewise small y 2 range such that Vs < Dat for 
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Hi — V2 £ [~ —0.25, 0.25]A. This means physically that only rigid structures are allowed 
and most of the base pair fluctuations are discarded as they overcome the Morse plateau 
thus yielding too large contributions to the classical action. Note that the number of turns 
added in Fig. [ltd) (with respect to Tw = 10) is equal to the number of turns subtracted in 
Fig. nib). Nevertheless an interesting, substantial asymmetry is found between the strongly 
overwound (Fig. [11(d)) and the strongly underwound (Fig. \W°)) conformation: the latter 
displays an energy scale much closer than the former to the equilibrium case (Fig. [lie)). 

Summing up, the double helix appears to be shaped by the balance between two com- 
peting needs: allowing for those fluctuational openings which are vital to biological func- 
tions meanwhile forbidding those fluctuations which may become too large and energeti- 
cally costly. General considerations based on the mesoscopic Hamiltonian in Eq. (pj) sug- 
gest that twists conformations with Tw ~ 10 emerge as the most suitable for the relaxed 
molecule. While the argument proposed so far cannot discern any appreciable difference 
between slightly underwound (Tw = 9) and slightly overwound (Tw =11) conformations, 
an energetic advantage is expected for stronger underwound (Tw < 9) molecules over the 
stronger overwound (Tw > 11) ones. After establishing that the stacking potential in 
Eq. Q captures the essentials of the interactions in the helicoidal geometry I proceed to 
develop a more quantitative analysis in terms of the path integral formalism [60] which also 
incorporates the environmental conditions due to temperature and solvent. 

3. Path Integral 

The double helix base pairs permanently fluctuate around their equilibrium positions 
represented by y n ~ in Eq. ([I]). Fluctuations get larger by increasing temperature and 
produce temporary openings, bubbles, which are peculiar of the biological functions [UJ H21 
HI]. Such fluctuational openings initiate at specific sites and propagate along the molecule 
backbone as a consequence of cooperative effects in the stacking. 

In an effort to account for the temperature dependent DNA dynamics, I have proposed 
[32] to apply path integral techniques to the mesoscopic Hamiltonian in Eq. Q. This 
permits to investigate the sequence specificity in the molecule and also to model the base 
pair fluctuational effects in the computation of the partition function. To accomplish a finite 
temperature description, the appropriate theoretical tool is the imaginary time path integral 



method [61] where the imaginary time r, an analytic continuation of the real time, runs on 
a scale which is set by the inverse temperature: r G [0, /3] with f3 = (/c^T) -1 . The length (3 
is partitioned in N T intervals each identified by a specific Tj. As this partitioning is repeated 
at any temperature while N T is kept constant, each Tj can vary inside the range whose upper 
bound is /3. Accordingly the Euclidean action integral transforms into a discrete sum (over 
i) of Hamiltonian contributions (evaluated at different r»). 

Thus, the space-time mapping technique for the model in Eq. (pi) assumes that the N 
base pair radial displacements y n are thought of as one dimensional paths x(tj) which are 
periodic functions of Tj, xfa) = x{ji + (3). % numbers the base pairs along the r-axis and 
Tj G [0, 0\, As a;(0) = x(/3), a molecule configuration is given by N T = N paths and, in 
the discrete imaginary time lattice, the separation between nearest neighbors base pairs is 
At = p/N T . 

Then, the displacements in Eq. ([I]) transform as: 



y n -> x{ji) 

y n _i ->x(n- Ar) (2) 



The base pair paths are written in Fourier series with cutoff Mp 



M F 

x(Ti) = x + ^2 \ a m cos(uj m Ti) + b m sin(c; m r i ' 



m=l 



u m = 2mir/f3 . (3) 

Accordingly, a set of coefficients {x , a m , b m } univocally defines a molecule configuration 
at a specific T, with the N base pairs given by {xfc), i = 1 , .., N T }. As a molecule state 
is defined by a point in the paths configuration space, path integration amounts to sum 
over a large ensemble of distinct molecule states each weighted by a Boltzmann factor. The 
method fully incorporates fluctuational effects around the ground state and monitors the 
bubble formation around and above room temperature. The classical partition function for 
the system in Eq. (PO) reads 
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Z c = <f>1)xexp[-PAc{x}] 

N T 

Ac{x} = Yl [^( r *) 2 + Vs(x(n),x(Ti - At)) + V M {x(n)) + V aol {x{n)) 

i=l 

(4) 
where the measure of integration T)x is expressed in terms of the Fourier coefficients by 

f 1 Mt * /m «2 /-At Mt 

iT)x = — =— / dx I K— x / da m / d& m . (5) 

7 V2A M J-ao, ^ v -V / ,/_ At j_ At 

Defining the free energy, F = j3~ 1 In Zc-, and the entropy, S = —dF/dT, the thermody- 
namical properties of the system can be computed via Eqs. Q, (J5J) taking a temperature 
range [T min ,T max ] in which fluctuational openings are expected to occur. The bounds, 
T min = 260K and T max = 390K, are chosen hereafter. 

Say N e ff the number of integration points in Eq. ^ to be selected on the base of 
thermodynamical and model potential physical requirements as detailed in Refs. [32j |33] . In 
particular, after setting an initial N e ff value at T m j n which ensures numerical convergence for 
Zq, the programme evaluates the free energy derivatives at any T > T min . If, for a given T, 
the entropy derivative is not positive then N e ff is increased and the thermodynamics is re- 
calculated. The procedure is reiterated for any T until the second law of thermodynamics is 
fulfilled throughout the whole temperature range. Eventually the method yields a N e ff(T) 
plot which denotes the number of possible trajectories for the i — th base pair whereas 
iV r x N e ff represents the total number of paths contributing to Zc- N T x N e ff yields the 
measure of the path ensemble size which, consistently with the entropy law, turns out to be 
a growing function of T. 



A M = ^tt/(3K, is the thermal wavelength in the classical regime. The cutoffs A^ and A^ 
in Eq. (fsl) remind us that the path amplitudes are temperature dependent. They are derived 
analytically noticing that the integration measure normalizes the free particle action |62j : 



Dx(r)exp - / rfr^x(r) 2 I. (6) 
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Hence, inserting Eqs. (J3]), ([5]), the l.h.s. of Eq. ^ transforms into a product of Gaufiian 
integrals whose normalization conditions yield: 

A T = Xfj,/ v 2 

mn 6 > 2 
Thus the temperature effects in the molecule states ensemble is included in the compu- 
tation also through the path amplitude cutoffs with oc vT dependence. The dimensionless 
U can be tuned to check the relevance of cutoff effects on the system thermodynamics [63J . 

4. Results 

The model is applied to a fragment with N T = 100 base pairs. The first Tj=i, ...,Ti=4s 
sites from left have the same sequence as the L48AS fragment of Refs.[6l] which shows 
fluctuations in the AT-rich side triggering the opening of distant GC pairs. In fact the 
melting curves as obtained by UV absorption measurements show that the base pairs open 
essentially in two stages, at two distinct temperatures. However, the fraction of base pairs 
opening in the first stage is larger than the whole fraction of AT- base pairs in the sequence. 
This implies that also GC- base pairs feel the effects of neighboring AT-rich regions, in 
particular of TATA boxes. 

These long range effects make the sequence of interest to us. Next, 52 AT-base pairs 
are added to L48AS on the right side to further develop fluctuational openings pS]- While 
the fragment has been studied in previous works [23 HZ] , it is here reconsidered in order to 
analyze the anharmonic effects in twisted geometries. 

The Tj=ioi site closes the double helix into a loop and therefore hosts a GC base pair. 
Inside AT-rich regions, pyrimidine-purine stacking steps favor local unwinding hence TA/TA 
steps are more efficient than AA/TT in driving bubble formation [4~3l |4"4"1 166] . The whole 
double strand sequence is: 

GC + 6AT + GC + 13AT + 8GC + AT + AGC + 
AT + AGC + AT + 8GC + [49 - 100]AT + GC . 

(8) 
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4.1 Melting Profiles 

As the UV signal changes quite abruptly when base pairs dissociate, the fraction of open 
base pairs is generally defined in terms of the Heaviside step function ■$(•) as: 

1 Nt 

f= jtY,*(< x <r*) > -Q > ( 9 ) 

T i=l 

where < x(t,j) >, the ensemble average for the displacement of the % — th base pair, is 
computed via Eq. Q with potential parameters as above. The threshold ( yields a criterion 
to establish whether an average base pair displacement is open, < x{t,j) >> (, or not. 

It should be remarked that Eq. (J9J) does not define univocally, on the T axis, at which 
stage the strand separation is large enough to produce full denaturation but, it provides a 
tool to evaluate quantitatively the amplitude of the fluctuational openings and the number 
of base pairs which participate in the bubble formation. While the choice of ( is somewhat 
arbitrary, in view of the fact that the B-DNA equilibrium diameter is ~ 20A, by taking 
( in the range ~ [0.5 — 2] A [67H69] one confidently explores the formation of nonlinear 
fluctuational openings in the appropriate temperature domain. 

The melting profiles given by Eq. (J9T) are plotted in Fig. [2] for moderately overtwisted 
and untwisted conformations: two turns of twist more (in Fig. [2|a)) and less (in Fig. [21(b)) 
respectively than (Tw) eq . Three ( values are chosen to calculate / as a function of T in the 
denaturation range. There is a remarkable effect due to the twist. For Tw = 8, all average 
displacements (/ = 1) get larger than C = 0.6A at T ~ 325K and ~ 30% (/ ~ 0.3) 
get larger than ( = lA in the upper T range. For Tw = 12, / ~ 1 for ( = 0.6A only at 
temperatures around 400f^ whereas, in such range, only ~ 15% of the average displacements 
exceed the ( = lA threshold. Overtwisting confers rigidity to the molecule while moderate 
untwisting provides flexibility. 

In both cases however segments of closed base pairs persist well inside the denaturation 
regime as indicated by recent neutron scattering measurements on B-DNA [70] . On the 
other hand, the Tw = 8 conformation shows that also the GC-base pairs (26% of the total) 
exceed the threshold ( = 0.6A at T ~ 325i^. In particular, the second GC-pair from left (in 
Eq. ([8])) has a large T dependent average displacement being embedded in a broad AT-rich 
region. This is consistent with the experimental findings on L48AS cited above. 
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FIG. 2: (Color online) Fractions of base pairs larger than 0.6, 0.8, 1A versus temperature in (a) 
moderately overwound, (b) moderately underwound double helix, (c) Total number of paths con- 
tributing to the partition function for four twist (and the ladder) conformations. The model 
potential parameters are as in Fig. [T] 



In Fig. 0(b), three main sharp enhancements appear at T ~ 280K, T ~ 325K ( £ = 
0.6 A) and T ~ 350^ ( ( = 0.8l ). Two smaller steps are at T ~ 300K ( ( = 0.6, I A 
plots). These are signatures of multistep openings around and above room temperature but 
~ 70% of the average displacements are still smaller than ( = lA at T ~ 400K suggesting 
that moderate untwisting is consistent with an overall thermal stability of the molecule. 

The Tw = 8 conformation is then expected to be similar, in terms of energetic conve- 
nience, to the (Tw) eq case. This is indeed confirmed by Fig. |2J(c) where N T x N e ff is plotted 
in the high T range. 

As explained in Section 3, iV T x N e ff is the number of possible paths which have been 
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selected to describe the base pair displacements. These paths are included in the compu- 
tation of Zq- While N T x N e ff varies with temperature, for a given T, it also depends on 
the twist geometry. The lower the number of paths required in the computation of the ther- 
modynamical properties, the higher is the stability of the twist conformation. N T x iV e // 
is in fact a measure of the system capability to sustain thermal fluctuations allowing for 
local openings without fully denaturating. The computation includes ensembles of order 10 7 
paths. The DPB ladder model and the strongly untwisted (Tw = 5) conformation display 
the largest N T x N e ff pointing to a major instability degree whereas relaxed (Tw = 10) 
and moderately untwisted (Tw = 8) conformations are the stablest ones. In fact the latter 
two plots overlap. Lying between the unstable and stable conformations, the overtwisted 
(Tw = 12) molecule shows a larger (than the Tw = 8, 10 cases) gradient in the upper T 
range. 

In summary, the analysis of the path ensemble size as a function of twist shows that a 
moderate unwinding of the double helix provides overall thermal stability to the molecule 
consistently with the computation of the average base pair displacements which determine 
the melting profiles. 

4.2 Anharmonic Stacking 

I investigate the interplay between twisting and stacking anharmonicity also in view of 
the special role given to the latter in the DPB model where a finite p induces a sharp 
denaturation transition driven by sizeable melting entropy (35]. Taking a ratio a/a n ~ 0.1, 
the range of the entropic barrier is here assumed to be longer than that of the Morse potential 
[T7] . In Fig. [3] the melting profiles of the DPB ladder model are calculated by increasing p 
over the previously used value, p = 2. Even very large p produce scant variations in the 
denaturation patterns pointing to a substantial irrelevance of the p driven anharmonicity for 
this conformation. Comparing Fig. |3Va) to Figs. W( a), (b), it is seen that the base pair opening 
in the ladder model is shifted at lower T, with / ~ 1 (for ( = 0.8A) already at T ~ 280K. 
This confirms that, in the absence of twist, the model is largely unstable consistently with 
the Tw = plot in Fig. |2J(c) and with the expectations provided by Fig. pTa). 

Quite different is the physical picture emerging from Fig. [4] where the efficacy of the 
stacking anharmonicity is sampled on the (Tw) eq conformation. Here, even slight enhance- 
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FIG. 3: (Color online) Fractions of average displacements larger than 0.6(O); 0.8(D), 1(§)A 
versus temperature in the DPB ladder model. Three anharmonic stacking p are assumed. The 
other stacking potential parameters are as in Fig. [TJ 

merits over the p = 2 value i) shift upwards along T the opening of the average base pair 
displacements and ii) flatten the melting profiles suggesting that the denaturation becomes 
more gradual. Note that p values in Figs. Efb),(c) are two orders of magnitude smaller than 
in Figs. [31(b), (c) respectively. The p = 8 plot says that, even at T ~ 400 K, there are no 
average displacements larger than lA while only 35% are larger than 0.6A Thus it is the 
twisting that switches on the p effect in the model. Anharmonic stacking induces those 
cooperative interactions along the molecule backbone which are peculiar of the ffuctuational 
openings [32] ■ At the same time anharmonic stacking renders the double helix flexible hence 
it does increase the molecule resilience against the whole thermal disruption of the hydrogen 
bonds. In this sense, anharmonicity is a stabilizing factor for the double helix. 



4.3 Solvent Potential 

Finally I discuss the effect of the solvent potential in Eq. ([I]) extending the range well 
beyond the length l s = 3 A used so far. l s defines the spatial range within which the sol- 
vent modifies the plateau of the Morse potential. A study of the DPB model [71] with 
the same potential term as in Eq. Q, concludes that growing l s would cause a sharp, first 
order-like, melting transition which is however smoothed by inclusion of a small twist angle. 
The melting profiles and the entropy plots shown in Fig. [5] do not support such conclu- 



16 




280 320 360 
T(K) 



280 320 360 
T(K) 



280 320 360 
T(K) 



FIG. 4: (Color online) As in Fig. KJ] but for the equilibrium twist conformation. In (b) and (c), p 
is smaller by a factor 100 than in Fig. p[b),(c) respectively. 



sions. Broadening the potential range shifts the denaturation steps upwards (Fig. [5ta)) but 
it does not change the smooth character of the denaturation transition which is moreover 
independent of the torsional degree: whatever twist angle is chosen the entropy grows con- 
tinuously. The curves for the equilibrium twist conformation are given in Fig. [5tb) both for 
l s = 3A and l s = 10A. In the latter case the entropy is slightly reduced consistently with 
the physical expectations attributing an overall stabilizing role to the solvent [72j . Instead, 
the entropy grows by enhancing the cutoff U in Eq. (JTl) due to the fact that much larger path 
amplitudes are included in the computation but even this does not sharpen the denaturation 
which remains continuous, second order-like. At the denaturation steps the melting entropy, 
~ 10 meVK~ l , is smaller by a factor 100 than that found in the DPB ladder model by 
transfer integral method [38] albeit for homogeneous DNA 



5. Discussion 

The path integral formalism has been applied to a circular heterogeneous DNA molecule 
modeled by a generalized Dauxois-Peyrard-Bishop Hamiltonian which incorporates a solvent 
potential and accounts for the helicoidal geometry. Inter-strand and intra-strand interac- 
tions are nonlinear. Sampling about 10 7 molecule states in the path configuration space, I 
have studied the system thermodynamics focussing on the interplay between stacking an- 
harmonicity and DNA supercoiling. The latter is tuned by varying the helical repeat with 
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FIG. 5: (Color online) (a) Fractions of average displacements larger than 0.6, 0.8, 1A in the equi- 
librium twist configuration for a large solvent potential range, l s = 10A The stacking potential 
parameters are as in Fig. p[a). (b) Equilibrium twist entropy as a function of i) l s and ii) the 
dimensionless path amplitude cutoff U. 

respect to the equilibrium conformation of room temperature B-DNA. Twisting has the ef- 
fect to bring closer to each other non-adjacent bases in the stack thus accounting for long 
range interactions which shape the molecule dynamics. Computations of melting profiles 
and size of the path ensemble consistently indicate that moderately untwisted helices, as 
much as the equilibrium conformation, are energetically favored with respect to the over- 
twisted ones. Then the path integral analysis quantitatively confirms the predictions based 
on the energetics of the stacking potential model versus degree of twist. 

At variance with previous studies, it is found that anharmonic stacking has little weight 
in the DPB model with zero twist whereas it becomes effective in helicoidal geometries and 
mainly for twist angles around the equilibrium conformation. In the latter an enhanced an- 
harmonicity smooths the melting profiles and shows an overall stabilizing role. Thus, I put 
forward a picture in which the nonlinear stacking potential is essential to make the molecule 
flexible and therefore able to sustain those thermal fructuational openings which are key to 
the biological functioning. This points to an interesting relation between intrinsic nonlinear 
character of the microscopic interactions and molecule topology. Accordingly anharmonicity 
is negligible in the ladder conformation which can be easily disordered by disruption of the 
hydrogen bonds. The calculations have been carried out, for a molecule with N = 100 base 
pairs, in a broad temperature range hosting several denaturation steps which start from the 
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AT-rich regions and spread to the GC-base pairs signalling the importance of cooperative 
effects in twisted geometries. While the melting patterns and the temperature location of 
the denaturation steps may change for other sequences with different relative weight of GC- 
and AT- base pairs, our findings regarding the interplay between anharmonicity and twist 
are independent of the specific sequence. Melting of heterogeneous DNA generally appears 
as an overall gradual crossover with presence of closed base pairs portions well inside the 
denaturation regime and with small entropy jumps. This is largely due to the strong fluctu- 
ational effects accounted for by the path integral approach while also heterogeneity smears 
the transition. The smooth nature of the denaturation is unaltered by the environmental 
conditions related to the solvent. 
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